function [output] = segmentationL(Image, modelHS)

numLabels = 4;

% Smooth/blur the image because of pixelation due to small image sizes
G = fspecial('disk',10);
Ifilt = imfilter(Image,G,'replicate');


% Use texture?
dummyMask = ones(size(Image,1),size(Image,2));
G = Image(:,:,2);
t1 = Run_MeasurementA(dummyMask,G);
t2 = Run_MeasurementB(dummyMask,G);

% Perform kmeans labeling using hue and saturation values
HSVfilt = rgb2hsv(Ifilt);
Hfilt = HSVfilt(:,:,1);
Sfilt = HSVfilt(:,:,2);
points = horzcat(Hfilt(:),Sfilt(:));
% k-means color segmentation
labeled = kmeansLabel(points,numLabels,'none');
labeled = reshape(labeled,size(dummyMask));


HSV = rgb2hsv(Image);
H = HSV(:,:,1);
S = HSV(:,:,2);

numCandidates = 0;
CCs = cell(numLabels,1);
for i=1:numLabels
    singleColor = (labeled==i);
    temp = bwareaopen(~singleColor, round(numel(singleColor)*0.005));
    singleColor = ~temp;
    cc = bwconncomp(singleColor);
    CCs{i} = cc;
    numCandidates = numCandidates + cc.NumObjects;
end

ROImasks = zeros(size(labeled,1),size(labeled,2),numCandidates);
numKeep = 0;
for i=1:numLabels
    cc = CCs{i};
    
    for j=1:cc.NumObjects
        candMask = zeros(size(labeled));
        candMask(cc.PixelIdxList{j}) = 1;
        hues = H(cc.PixelIdxList{j});
        sats = S(cc.PixelIdxList{j});
        
        s = regionprops(cc,'Area');
        area = cat(1, s.Area);
        area = area(j);

        % only consider regions greater than certain size
        if area/(size(candMask,1)*size(candMask,2)) > 0.005;
            if mean(hues) <= mean(sats)
                sumHS = sum([mean(hues) mean(sats)]);
                dists = abs(sum(modelHS,2)-sumHS);
                if min(dists) < std(sum(modelHS,2))
                    numKeep = numKeep + 1;
                    ROImasks(:,:,numKeep) = candMask;
                end
            end
        end
    end  
end


if numKeep > 0
    %figure,imshow(labeled,[])
%     subplot(1,numKeep,1)
%     for k=1:numKeep
%         subplot(1,numKeep,k)
%         imshow(ROImasks(:,:,k));
%     end
%     clf
    output = ROImasks(:,:,1:numKeep);
else
    %Labels of corner pixels provide information about potential background
    %regions
    G = fspecial('disk',10);
    Ifilt = imfilter(Image,G,'replicate');
    imshow(Ifilt);
    
    labeled = kmeansColor(Ifilt,4,'rgb');
    
    [m n] = size(labeled);
    corners = [labeled(1,1), labeled(1,n), labeled(m,1), labeled(m,n)];
    mo1 = mode(corners);
    corners(corners==mo1) = [];
    mo2 = mode(corners);
    
    mask = labeled~=mo1 & labeled~=mo2;
    mask = imfill(mask,'holes');
    output = mask;
end





% Labels of corner pixels provide information about potential background
% regions
% [m n] = size(labeled);
% corners = [labeled(1,1), labeled(1,n), labeled(m,1), labeled(m,n)];
% mo1 = mode(corners); 
% corners(corners==mo1) = []; 
% mo2 = mode(corners);
% 
% mask = labeled~=mo1 & labeled~=mo2;
% mask = imfill(mask,'holes');

% G = fspecial('disk',10);
% Ifilt = imfilter(Image,G,'replicate');
% imshow(Ifilt);
% 
% labeled = kmeansColor(Ifilt,5,'rgb');
% imshow(labeled, []);
% 
% segmented_images = cell(1,6);
% 
% count = 0;
% for k = 1:6
%     color = labeled;
%     color(labeled ~= k) = 0;
%     segmented_images{k} = color;
%     segmented_images{k} = imclearborder(segmented_images{k});
%     if sum(segmented_images{k}(:)) > 0
%         count = count+1;
%     end
% end
% 
% mask = segmented_images{1}|segmented_images{2}|segmented_images{3}|segmented_images{4}|segmented_images{5}|segmented_images{6};
% if count < 3
%     labeled = kmeansColor(Ifilt,3,'rgb');
%     disp('redo k means');
%     [m n] = size(labeled);
%     corners = [labeled(1,1), labeled(1,n), labeled(m,1), labeled(m,n)];
%     mo1 = mode(corners); corners(corners==mo1) = []; mo2 = mode(corners);
%     if ~isnan(mo2)
%         mask = labeled~=mo1 & labeled~=mo2;
%     else
%         mask = labeled~=mo1;
%     end
%     mask = imfill(mask,'holes');
%     output(:,:,1) = mask;
% else
%     mask = medfilt2(mask,[6 6]);
%     se = strel('disk',4);
%     mask = imfill(mask,'holes');
%     mask = imdilate(mask, se);
%     mask = imfill(mask,'holes');
%     mask = imerode(mask, se);
% 
%     cc = bwconncomp(mask);
% 
%     [row col] = size(mask);
% 
%     s = regionprops(cc,'Area');
%     area = cat(1, s.Area);
% 
%     keep = area/(size(mask,1)*size(mask,2)) > 0.0001;
%     newNumObjects = sum(keep);
% 
%     for i=1:cc.NumObjects
%         if ~keep(i)
%             mask(cc.PixelIdxList{i}) = 0;
%         end
%     end
% 
%     output = zeros(row, col, newNumObjects);
%     output = logical(output);
% 
%     if newNumObjects > 1
%         for i = 1:cc.NumObjects
%             if keep(i)
%                 label = zeros(size(mask));
%                 label(cc.PixelIdxList{i}) = 1;
%                 output(:,:,i) = label;
%             end
%         end
%     else
%         output(:,:,1) = mask;
%     end
% end